b.vec <-  c(.25, .5, .75, 1, 1.25)[S]  # Coefficient of explanatory variable
s <- 1000 # Number of Simulations
result.array <- array(NA, c(6,8,s))
ptm <- proc.time()
for(i in 1:s){
print(paste("Simulation #", S, i, sep=""))
#######################################################################################################################
### Simulate Data Set
#c.vec <-  rep(0, k) # Random correlations
#M <- huerlimann(k, c.vec)
M <- diag(k)
X <- rmvnorm(n, mean=rep(0, k), sigma = M)
X### DGP
Ey <- X[, 1] * b.vec
res <- rnorm(n, 0, 8) # Fix this bit to get SE of around .25 !
y <- Ey + res
fake.data <- as.data.frame(cbind(y, X))
colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10")
#colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10",
#                        "x11", "x12", "x13", "x14", "x15")
#                         ,"x16", "x17", "x18", "x19", "x20")
#######################################################################################################################
#### Generate Variable Combinations
X <- fake.data[,-1]
combs <- do.call(expand.grid, rep(list(c(FALSE, TRUE)), ncol(X)))[-1,]
### Assign Variable Names
vars <- apply(combs, 1, function(i) names(X)[i])
form <- paste("y ~ ", lapply(vars, paste, collapse = " + "), sep = "")
form <- lapply(form, as.formula)
### Fit Models
#ptm <- proc.time()
models <- lapply(form, fitmodel) # fit models
#proc.time() - ptm
### Post-Process Models
names(models) <- seq_along(models) # assign model ID
mods <- mod_summary(models) # summary of models
coefs <- coef_summary(models) # summary of coefficients
coefs_wide_b <- acast(coefs, model ~ variable, value.var=c("Est"))
coefs_wide_se <- acast(coefs, model ~ variable, value.var=c("SE"))
colnames(coefs_wide_se) <- paste("SE", colnames(coefs_wide_b) , sep="")
coefs <- cbind(coefs_wide_b, coefs_wide_se)
coefs <- coefs[,order(colnames(coefs), decreasing=T)]
model <- c(1:dim(mods)[1])
coefs <- cbind(model, coefs)
DATA <- merge(mods, coefs, by="model")
#DATA <- DATA[DATA[,"df"]==2, ]
###############################################################################
# EBA vs. Model Averaging
K <- k
Q <- q
true.b <- c(b.vec, rep(0, K-Q))
## EBA (Leamer)
eb.min <- apply(DATA[,7:(6+K)], 2, min, na.rm=T)
eb.max <- apply(DATA[,7:(6+K)], 2, max, na.rm=T)
index.min <- apply(DATA[,7:(6+K)], 2, function(i){which(i==min(i,na.rm=T))})
index.max <- apply(DATA[,7:(6+K)], 2, function(i){which(i==max(i,na.rm=T))})
if(class(index.min)=="list"){
index.min <- unlist(lapply(index.min, min))
}
if(class(index.max)=="list"){
index.max <- unlist(lapply(index.max, max))
}
se.min <- DATA[,(7+K):(7+K+K-1)][cbind(index.min, c(1:K))]
se.max <- DATA[,(7+K):(7+K+K-1)][cbind(index.max, c(1:K))]
# Dreher
qu.1 <- apply(DATA[,7:(6+K)], 2, quantile, .1, na.rm=T)
qu.9 <- apply(DATA[,7:(6+K)], 2, quantile, .9, na.rm=T)
# Model Averaging
ma.ew <- apply(DATA[,7:(6+K)], 2, mean, na.rm=T)
var.b.ew <- apply(DATA[,(7+K):(7+K+K-1)], 2, mean, na.rm=T)
var.b.ew <- var.b.ew^2
# var.ew <- apply(DATA[,7:(6+K)], 2, var, na.rm=T)
var.ew <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 , na.rm=T)/sum(i, na.rm=T)  })
sd.ew <- sqrt(var.b.ew + var.ew)
bic <- DATA$BIC/1000 # rescale BIC
w <- (exp(-bic/2))/(sum(exp(-bic/2))) # weight (cf. Moral-Benito 2010: 19)
ma.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.bic <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
var.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 * w, na.rm=T)})
sd.w.bic <- sqrt(var.b.w.bic + var.w.bic)
dev <- DATA$Deviance/1000 # rescale Deviance
L <- exp(-dev/2) # Likelihood
w <- L/sum(L) # weight
ma.w.ll <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.ll <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
sd.w.ll <- sqrt(var.b.w.ll)
############################################################################################################################
# Robustness Criteria
in.out.eb.s <- 0<eb.min*eb.max & 1.96<abs(eb.min/se.min) & 1.96<abs(eb.max/se.max) # EBA: min and max same sign & significant
in.out.eb <- 0<eb.min*eb.max # EBA: min and max same sign
in.out.ma <- 1.96<abs(ma.ew/sd.ew) # Unweighted Model Averaging
in.out.ma.sim <- .95<apply(cbind(pnorm(0, ma.w.ll, sd.w.ll), 1-pnorm(0, ma.w.ll, sd.w.ll)), 1, max) # Sala-i-Martin: Model Averaging weighted by LL
in.out.ma.bic <- 1.96<abs(ma.w.bic/sd.w.bic) # Model Averaging weighted by BIC
in.out.d90 <- 0<qu.1*qu.9 # Dreher: 90% quantiles same sign
beta.mat <- rbind(in.out.eb.s, in.out.eb, in.out.ma, in.out.ma.sim, in.out.ma.bic, in.out.d90, ma.ew, ma.w.ll, ma.w.bic, qu.1, eb.min)
colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 10, 1)
#colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 14, 13, 12, 11, 10, 1)
beta.mat <- beta.mat[,order(as.numeric(colnames(beta.mat)), decreasing=F)]
#############################################################################################################################
# Predictive Criteria
num.robust <- apply(beta.mat[1:6, 1:K], 1, sum) # Number of variables declared "robust"
true.pos <- beta.mat[1:6, 1:Q] # Number of True Positives
false.pos <- apply(beta.mat[1:6, (Q+1):K], 1, sum) # Number of FP
true.neg <- (K-Q)-false.pos
false.neg <- Q-true.pos
tpr <- true.pos/(true.pos+false.neg) # TPR
tnr <- true.neg/(true.neg+false.pos) # TNR
correct.sign.ma.ew <- sum(beta.mat["ma.ew",1:Q]>0)/sum(true.b[1:Q]>0)  ### Correct Sign
correct.sign.ma.w.ll <- sum(beta.mat["ma.w.ll",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.ma.w.bic <- sum(beta.mat["ma.w.bic",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.d90 <-  sum(beta.mat["in.out.d90", 1:Q]==1 & beta.mat["qu.1", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb <-  sum(beta.mat["in.out.eb", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb.s <-  sum(beta.mat["in.out.eb.s", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign <- c(correct.sign.eb.s, correct.sign.eb, correct.sign.ma.ew, correct.sign.ma.w.ll, correct.sign.ma.w.bic, correct.sign.d90)
rmse.ew <- sqrt(mean((true.b - beta.mat["ma.ew",])^2)) ## RMSE
rmse.w.ll <- sqrt(mean((true.b - beta.mat["ma.w.ll",])^2))
rmse.w.bic <- sqrt(mean((true.b - beta.mat["ma.w.bic",])^2))
rmse <- c(NA, NA, rmse.ew, rmse.w.ll, rmse.w.bic, NA)
corr.ew <- cor(true.b, beta.mat["ma.ew",]) ## Correlation with true b
corr.w.ll <- cor(true.b, beta.mat["ma.w.ll",])
corr.w.bic <- cor(true.b, beta.mat["ma.w.bic",])
corrs <- c(NA, NA, corr.ew, corr.w.ll, corr.w.bic, NA)
result.mat <- cbind(num.robust, true.pos, false.pos, true.neg, false.neg, correct.sign, rmse, corrs)
result.array[,,i] <- result.mat
}
time.needed <- proc.time()-ptm
colnames(result.array) <- c("num.robust", "true.pos", "false.pos", "true.neg", "false.neg", "correct.sign", "rmse", "corrs")
rownames(result.array) <- c("eb.s", "eb", "ma", "ma.sim", "ma.bic", "d90")
saveRDS(result.array, paste("one_result_mat_", S, ".rds", sep=""))
}
for(S in 1:5){
n <- 1000
k <- 10
q <- 1
r2 <- .5
b.vec <-  c(.25, .5, .75, 1, 1.25)[S]  # Coefficient of explanatory variable
s <- 1000 # Number of Simulations
result.array <- array(NA, c(6,8,s))
ptm <- proc.time()
for(i in 1:s){
print(paste("Condition: ", S, " Simulation: ", i, sep=""))
#######################################################################################################################
### Simulate Data Set
#c.vec <-  rep(0, k) # Random correlations
#M <- huerlimann(k, c.vec)
M <- diag(k)
X <- rmvnorm(n, mean=rep(0, k), sigma = M)
X### DGP
Ey <- X[, 1] * b.vec
res <- rnorm(n, 0, 8) # Fix this bit to get SE of around .25 !
y <- Ey + res
fake.data <- as.data.frame(cbind(y, X))
colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10")
#colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10",
#                        "x11", "x12", "x13", "x14", "x15")
#                         ,"x16", "x17", "x18", "x19", "x20")
#######################################################################################################################
#### Generate Variable Combinations
X <- fake.data[,-1]
combs <- do.call(expand.grid, rep(list(c(FALSE, TRUE)), ncol(X)))[-1,]
### Assign Variable Names
vars <- apply(combs, 1, function(i) names(X)[i])
form <- paste("y ~ ", lapply(vars, paste, collapse = " + "), sep = "")
form <- lapply(form, as.formula)
### Fit Models
#ptm <- proc.time()
models <- lapply(form, fitmodel) # fit models
#proc.time() - ptm
### Post-Process Models
names(models) <- seq_along(models) # assign model ID
mods <- mod_summary(models) # summary of models
coefs <- coef_summary(models) # summary of coefficients
coefs_wide_b <- acast(coefs, model ~ variable, value.var=c("Est"))
coefs_wide_se <- acast(coefs, model ~ variable, value.var=c("SE"))
colnames(coefs_wide_se) <- paste("SE", colnames(coefs_wide_b) , sep="")
coefs <- cbind(coefs_wide_b, coefs_wide_se)
coefs <- coefs[,order(colnames(coefs), decreasing=T)]
model <- c(1:dim(mods)[1])
coefs <- cbind(model, coefs)
DATA <- merge(mods, coefs, by="model")
#DATA <- DATA[DATA[,"df"]==2, ]
###############################################################################
# EBA vs. Model Averaging
K <- k
Q <- q
true.b <- c(b.vec, rep(0, K-Q))
## EBA (Leamer)
eb.min <- apply(DATA[,7:(6+K)], 2, min, na.rm=T)
eb.max <- apply(DATA[,7:(6+K)], 2, max, na.rm=T)
index.min <- apply(DATA[,7:(6+K)], 2, function(i){which(i==min(i,na.rm=T))})
index.max <- apply(DATA[,7:(6+K)], 2, function(i){which(i==max(i,na.rm=T))})
if(class(index.min)=="list"){
index.min <- unlist(lapply(index.min, min))
}
if(class(index.max)=="list"){
index.max <- unlist(lapply(index.max, max))
}
se.min <- DATA[,(7+K):(7+K+K-1)][cbind(index.min, c(1:K))]
se.max <- DATA[,(7+K):(7+K+K-1)][cbind(index.max, c(1:K))]
# Dreher
qu.1 <- apply(DATA[,7:(6+K)], 2, quantile, .1, na.rm=T)
qu.9 <- apply(DATA[,7:(6+K)], 2, quantile, .9, na.rm=T)
# Model Averaging
ma.ew <- apply(DATA[,7:(6+K)], 2, mean, na.rm=T)
var.b.ew <- apply(DATA[,(7+K):(7+K+K-1)], 2, mean, na.rm=T)
var.b.ew <- var.b.ew^2
# var.ew <- apply(DATA[,7:(6+K)], 2, var, na.rm=T)
var.ew <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 , na.rm=T)/sum(i, na.rm=T)  })
sd.ew <- sqrt(var.b.ew + var.ew)
bic <- DATA$BIC/1000 # rescale BIC
w <- (exp(-bic/2))/(sum(exp(-bic/2))) # weight (cf. Moral-Benito 2010: 19)
ma.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.bic <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
var.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 * w, na.rm=T)})
sd.w.bic <- sqrt(var.b.w.bic + var.w.bic)
dev <- DATA$Deviance/1000 # rescale Deviance
L <- exp(-dev/2) # Likelihood
w <- L/sum(L) # weight
ma.w.ll <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.ll <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
sd.w.ll <- sqrt(var.b.w.ll)
############################################################################################################################
# Robustness Criteria
in.out.eb.s <- 0<eb.min*eb.max & 1.96<abs(eb.min/se.min) & 1.96<abs(eb.max/se.max) # EBA: min and max same sign & significant
in.out.eb <- 0<eb.min*eb.max # EBA: min and max same sign
in.out.ma <- 1.96<abs(ma.ew/sd.ew) # Unweighted Model Averaging
in.out.ma.sim <- .95<apply(cbind(pnorm(0, ma.w.ll, sd.w.ll), 1-pnorm(0, ma.w.ll, sd.w.ll)), 1, max) # Sala-i-Martin: Model Averaging weighted by LL
in.out.ma.bic <- 1.96<abs(ma.w.bic/sd.w.bic) # Model Averaging weighted by BIC
in.out.d90 <- 0<qu.1*qu.9 # Dreher: 90% quantiles same sign
beta.mat <- rbind(in.out.eb.s, in.out.eb, in.out.ma, in.out.ma.sim, in.out.ma.bic, in.out.d90, ma.ew, ma.w.ll, ma.w.bic, qu.1, eb.min)
colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 10, 1)
#colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 14, 13, 12, 11, 10, 1)
beta.mat <- beta.mat[,order(as.numeric(colnames(beta.mat)), decreasing=F)]
#############################################################################################################################
# Predictive Criteria
num.robust <- apply(beta.mat[1:6, 1:K], 1, sum) # Number of variables declared "robust"
true.pos <- beta.mat[1:6, 1:Q] # Number of True Positives
false.pos <- apply(beta.mat[1:6, (Q+1):K], 1, sum) # Number of FP
true.neg <- (K-Q)-false.pos
false.neg <- Q-true.pos
tpr <- true.pos/(true.pos+false.neg) # TPR
tnr <- true.neg/(true.neg+false.pos) # TNR
correct.sign.ma.ew <- sum(beta.mat["ma.ew",1:Q]>0)/sum(true.b[1:Q]>0)  ### Correct Sign
correct.sign.ma.w.ll <- sum(beta.mat["ma.w.ll",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.ma.w.bic <- sum(beta.mat["ma.w.bic",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.d90 <-  sum(beta.mat["in.out.d90", 1:Q]==1 & beta.mat["qu.1", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb <-  sum(beta.mat["in.out.eb", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb.s <-  sum(beta.mat["in.out.eb.s", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign <- c(correct.sign.eb.s, correct.sign.eb, correct.sign.ma.ew, correct.sign.ma.w.ll, correct.sign.ma.w.bic, correct.sign.d90)
rmse.ew <- sqrt(mean((true.b - beta.mat["ma.ew",])^2)) ## RMSE
rmse.w.ll <- sqrt(mean((true.b - beta.mat["ma.w.ll",])^2))
rmse.w.bic <- sqrt(mean((true.b - beta.mat["ma.w.bic",])^2))
rmse <- c(NA, NA, rmse.ew, rmse.w.ll, rmse.w.bic, NA)
corr.ew <- cor(true.b, beta.mat["ma.ew",]) ## Correlation with true b
corr.w.ll <- cor(true.b, beta.mat["ma.w.ll",])
corr.w.bic <- cor(true.b, beta.mat["ma.w.bic",])
corrs <- c(NA, NA, corr.ew, corr.w.ll, corr.w.bic, NA)
result.mat <- cbind(num.robust, true.pos, false.pos, true.neg, false.neg, correct.sign, rmse, corrs)
result.array[,,i] <- result.mat
}
time.needed <- proc.time()-ptm
colnames(result.array) <- c("num.robust", "true.pos", "false.pos", "true.neg", "false.neg", "correct.sign", "rmse", "corrs")
rownames(result.array) <- c("eb.s", "eb", "ma", "ma.sim", "ma.bic", "d90")
saveRDS(result.array, paste("one_result_mat_", S, ".rds", sep=""))
}
for(i in 1:s){
print(paste("Experiment: 1A", " Condition: ", S, " Simulation: ", i, sep=""))
#######################################################################################################################
### Simulate Data Set
#c.vec <-  rep(0, k) # Random correlations
#M <- huerlimann(k, c.vec)
M <- diag(k)
X <- rmvnorm(n, mean=rep(0, k), sigma = M)
X### DGP
Ey <- X[, 1] * b.vec
res <- rnorm(n, 0, 8) # Fix this bit to get SE of around .25 !
y <- Ey + res
fake.data <- as.data.frame(cbind(y, X))
colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10")
#colnames(fake.data) <- c("y", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10",
#                        "x11", "x12", "x13", "x14", "x15")
#                         ,"x16", "x17", "x18", "x19", "x20")
#######################################################################################################################
#### Generate Variable Combinations
X <- fake.data[,-1]
combs <- do.call(expand.grid, rep(list(c(FALSE, TRUE)), ncol(X)))[-1,]
### Assign Variable Names
vars <- apply(combs, 1, function(i) names(X)[i])
form <- paste("y ~ ", lapply(vars, paste, collapse = " + "), sep = "")
form <- lapply(form, as.formula)
### Fit Models
#ptm <- proc.time()
models <- lapply(form, fitmodel) # fit models
#proc.time() - ptm
### Post-Process Models
names(models) <- seq_along(models) # assign model ID
mods <- mod_summary(models) # summary of models
coefs <- coef_summary(models) # summary of coefficients
coefs_wide_b <- acast(coefs, model ~ variable, value.var=c("Est"))
coefs_wide_se <- acast(coefs, model ~ variable, value.var=c("SE"))
colnames(coefs_wide_se) <- paste("SE", colnames(coefs_wide_b) , sep="")
coefs <- cbind(coefs_wide_b, coefs_wide_se)
coefs <- coefs[,order(colnames(coefs), decreasing=T)]
model <- c(1:dim(mods)[1])
coefs <- cbind(model, coefs)
DATA <- merge(mods, coefs, by="model")
#DATA <- DATA[DATA[,"df"]==2, ]
###############################################################################
# EBA vs. Model Averaging
K <- k
Q <- q
true.b <- c(b.vec, rep(0, K-Q))
## EBA (Leamer)
eb.min <- apply(DATA[,7:(6+K)], 2, min, na.rm=T)
eb.max <- apply(DATA[,7:(6+K)], 2, max, na.rm=T)
index.min <- apply(DATA[,7:(6+K)], 2, function(i){which(i==min(i,na.rm=T))})
index.max <- apply(DATA[,7:(6+K)], 2, function(i){which(i==max(i,na.rm=T))})
if(class(index.min)=="list"){
index.min <- unlist(lapply(index.min, min))
}
if(class(index.max)=="list"){
index.max <- unlist(lapply(index.max, max))
}
se.min <- DATA[,(7+K):(7+K+K-1)][cbind(index.min, c(1:K))]
se.max <- DATA[,(7+K):(7+K+K-1)][cbind(index.max, c(1:K))]
# Dreher
qu.1 <- apply(DATA[,7:(6+K)], 2, quantile, .1, na.rm=T)
qu.9 <- apply(DATA[,7:(6+K)], 2, quantile, .9, na.rm=T)
# Model Averaging
ma.ew <- apply(DATA[,7:(6+K)], 2, mean, na.rm=T)
var.b.ew <- apply(DATA[,(7+K):(7+K+K-1)], 2, mean, na.rm=T)
var.b.ew <- var.b.ew^2
# var.ew <- apply(DATA[,7:(6+K)], 2, var, na.rm=T)
var.ew <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 , na.rm=T)/sum(i, na.rm=T)  })
sd.ew <- sqrt(var.b.ew + var.ew)
bic <- DATA$BIC/1000 # rescale BIC
w <- (exp(-bic/2))/(sum(exp(-bic/2))) # weight (cf. Moral-Benito 2010: 19)
ma.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.bic <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
var.w.bic <- apply(DATA[,7:(6+K)], 2, function(i){sum( (i-mean(i, na.rm=T))^2 * w, na.rm=T)})
sd.w.bic <- sqrt(var.b.w.bic + var.w.bic)
dev <- DATA$Deviance/1000 # rescale Deviance
L <- exp(-dev/2) # Likelihood
w <- L/sum(L) # weight
ma.w.ll <- apply(DATA[,7:(6+K)], 2, function(i){sum(i * w, na.rm=T)})
var.b.w.ll <- apply(DATA[,(7+K):(7+K+K-1)], 2, function(i){sum((i)^2 * w, na.rm=T)})
sd.w.ll <- sqrt(var.b.w.ll)
############################################################################################################################
# Robustness Criteria
in.out.eb.s <- 0<eb.min*eb.max & 1.96<abs(eb.min/se.min) & 1.96<abs(eb.max/se.max) # EBA: min and max same sign & significant
in.out.eb <- 0<eb.min*eb.max # EBA: min and max same sign
in.out.ma <- 1.96<abs(ma.ew/sd.ew) # Unweighted Model Averaging
in.out.ma.sim <- .95<apply(cbind(pnorm(0, ma.w.ll, sd.w.ll), 1-pnorm(0, ma.w.ll, sd.w.ll)), 1, max) # Sala-i-Martin: Model Averaging weighted by LL
in.out.ma.bic <- 1.96<abs(ma.w.bic/sd.w.bic) # Model Averaging weighted by BIC
in.out.d90 <- 0<qu.1*qu.9 # Dreher: 90% quantiles same sign
beta.mat <- rbind(in.out.eb.s, in.out.eb, in.out.ma, in.out.ma.sim, in.out.ma.bic, in.out.d90, ma.ew, ma.w.ll, ma.w.bic, qu.1, eb.min)
colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 10, 1)
#colnames(beta.mat) <- c(9, 8, 7, 6, 5, 4, 3, 2, 14, 13, 12, 11, 10, 1)
beta.mat <- beta.mat[,order(as.numeric(colnames(beta.mat)), decreasing=F)]
#############################################################################################################################
# Predictive Criteria
num.robust <- apply(beta.mat[1:6, 1:K], 1, sum) # Number of variables declared "robust"
true.pos <- beta.mat[1:6, 1:Q] # Number of True Positives
false.pos <- apply(beta.mat[1:6, (Q+1):K], 1, sum) # Number of FP
true.neg <- (K-Q)-false.pos
false.neg <- Q-true.pos
tpr <- true.pos/(true.pos+false.neg) # TPR
tnr <- true.neg/(true.neg+false.pos) # TNR
correct.sign.ma.ew <- sum(beta.mat["ma.ew",1:Q]>0)/sum(true.b[1:Q]>0)  ### Correct Sign
correct.sign.ma.w.ll <- sum(beta.mat["ma.w.ll",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.ma.w.bic <- sum(beta.mat["ma.w.bic",1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.d90 <-  sum(beta.mat["in.out.d90", 1:Q]==1 & beta.mat["qu.1", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb <-  sum(beta.mat["in.out.eb", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign.eb.s <-  sum(beta.mat["in.out.eb.s", 1:Q]==1 & beta.mat["eb.min", 1:Q]>0)/sum(true.b[1:Q]>0)
correct.sign <- c(correct.sign.eb.s, correct.sign.eb, correct.sign.ma.ew, correct.sign.ma.w.ll, correct.sign.ma.w.bic, correct.sign.d90)
rmse.ew <- sqrt(mean((true.b - beta.mat["ma.ew",])^2)) ## RMSE
rmse.w.ll <- sqrt(mean((true.b - beta.mat["ma.w.ll",])^2))
rmse.w.bic <- sqrt(mean((true.b - beta.mat["ma.w.bic",])^2))
rmse <- c(NA, NA, rmse.ew, rmse.w.ll, rmse.w.bic, NA)
corr.ew <- cor(true.b, beta.mat["ma.ew",]) ## Correlation with true b
corr.w.ll <- cor(true.b, beta.mat["ma.w.ll",])
corr.w.bic <- cor(true.b, beta.mat["ma.w.bic",])
corrs <- c(NA, NA, corr.ew, corr.w.ll, corr.w.bic, NA)
result.mat <- cbind(num.robust, true.pos, false.pos, true.neg, false.neg, correct.sign, rmse, corrs)
result.array[,,i] <- result.mat
}
source("PSRM_one_mcs.R")  # Runs Part 1 of the Simulations (p=1)
source("PSRM_five_runif_mcs.R") # Runs Part B of the Simulations (p=5)
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
source("PSRM_one_runif_mcs.R") # Runs Part A of the Simulations (p=1)
source("PSRM_five_runif_mcs.R") # Runs Part B of the Simulations (p=5)
source("PSRM_nine_runif_mcs.R") # Runs Part C of the Simulations (p=9)
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
source("PSRM_five_runif_mcs.R") # Runs Part B of the Simulations (p=5)
# Specify Working Directory
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
# Load Required R Packages
library(mvtnorm)
library(reshape2)
library(plyr)
# Load Custom Functions
source("huerlimann_function.R") # Function to simulate NxN correlation matrizes
source("extract_fun_OLS.R") # Function to extract model quantities
# Specify Working Directory
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
# Load Required R Packages
library(mvtnorm)
library(reshape2)
library(plyr)
# Load Custom Functions
source("huerlimann_function.R") # Function to simulate NxN correlation matrizes
source("extract_fun_OLS.R") # Function to extract model quantities
source("PSRM_one_mcs.R")  # Runs Part A of the Simulations (p=1)
source("PSRM_five_mcs.R") # Runs Part B of the Simulations (p=5)
source("PSRM_nine_mcs.R") # Runs Part C of the Simulations (p=9)
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
MeanOutFun <- function(IN_FILE, OUT_FILE){
result.array <- readRDS(paste("/Users/richardtraunmuller/Dropbox/rep_hegsamb/", IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
for(j in 1:8){
mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T)
}
}
write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
#print(mean.array)
}
p <- 1
n1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)
p1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))
tab.1.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])
n1
n1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n1
MeanOutFun <- function(IN_FILE, OUT_FILE){
result.array <- readRDS(paste("/Users/richardtraunmuller/Dropbox/rep_hegsamb/", IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
for(j in 1:8){
mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T)
}
}
write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
print(mean.array)
}
p <- 1
n1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)
p1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))
tab.1.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
source("PSRM_summary_2.R") # Produces a .csv File of Table 2 in the Paper
source("PSRM_summary_3.R") # Produces a .csv File of Table 3 in the Paper
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
source("PSRM_summary_2.R") # Produces a .csv File of Table 2 in the Paper
source("PSRM_summary_3.R") # Produces a .csv File of Table 3 in the Paper
setwd(paste(YOURPATH, "PSRM_Replication/Simulation_Results/", sep=""))
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
YOURPATH <- "/Users/richardtraunmuller/Dropbox/" # Change Your Path Here
setwd(paste(YOURPATH, "PSRM_Replication/", sep=""))
source("PSRM_summary_2.R") # Produces a .csv File of Table 2 in the Paper
source("PSRM_summary_1.R") # Produces a .csv File of Table 1 in the Paper
source("PSRM_summary_3.R") # Produces a .csv File of Table 3 in the Paper
source("PSRM_summary_2.R") # Produces a .csv File of Table 2 in the Paper
20/99237
3879/99237
95338/99237
737575+99237
99237/(737575+99237)
library(knitr)
knitr::opts_chunk$set(echo = FALSE)
mean(data$Attack)
mean(data$Attack, na.rm=T)
data <- readRDS("/Users/richardtraunmuller/Dropbox/footballdata.rds")
mean(data$Attack, na.rm=T)
